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Abstract 

Any nonspherical distribution of density inside planets and stars gives rise to a non-spherical 
external gravity and change of shape. If part or all of the observed zonal flows at the cloud 
deck of Jupiter and Saturn represent deep interior dynamics, then the density perturbations 
associated with the deep zonal flows could generate gravitational signals detectable by the 
Juno mission and the Cassini Grand Finale. Here we present a critical examination of the 
applicability of the thermal wind equation to calculate the wind induced gravity moments. 
Our analysis shows that wind induced gravity moments calculated from TWE are in overall 
agreement with the full solution to the Euler equation. However, the accuracy of individual 
high-degree moments calculated from TWE depends crucially on retaining the nonspheric¬ 
ity of the background density and gravity. Only when the background nonsphericity of the 
planet is taken into account, does the TWE make accurate enough prediction (with a few tens 
of percent errors) for individual high-degree gravity moments associated with deep zonal 
flows. Since the TWE is derived from the curl of the Euler equation and is a local relation, it 
necessarily says nothing about any density perturbations that contribute irrotational terms to 
the Euler equation and that have a non-local origin. However, the predicted corrections from 
these density contributions to the low harmonic degree gravity moments are not discernible 
from insignificant changes in interior models while the corrections at high harmonic degree 
are very small, a few percent or less. 

1 Introduction 

The interior structures and dynamics of the solar system giant planets remain elusive 
after decades of observational, experimental and theoretical studies [cf Stevenson, 1982; 
Hubbard et ai, 2002; Guillot, 2005; Guillot and Gautier, 2014, and reference therein]. Eor 
example, whether present-day Jupiter and Saturn have well-defined cores remains an open 
question; the total enrichment of heavy elements inside Jupiter and Saturn are not well con¬ 
strained; the structural and dynamical consequences of the likely on-going sedimentation of 
helium and neon inside Jupiter and Saturn have not been fully worked out [cf Stevenson and 
Salpeter, 1977; Fortney and Hubbard, 2003; Nettelmann et al., 2015]. 

One lasting debate concerns the nature of the observed east-west zonal flows on the 
cloud layers of giant planets with amplitude on the order of 100 mis: no consensus has been 
reached upon whether these zonal winds represent shallow atmospheric dynamics or deep 
interior dynamics [e.g. Vasavada and Showman, 2005; Liu et ai, 2008; Jones and Kuzanyan, 
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2009; Kaspi et al, 2009; Liu and Schneider, 2010; Gastine et ai, 2013]. The forward fluid 
dynamics and magnetohydrodynamics (MHD) problem about the nature of giant planet zonal 
flows may be hard to settle given the complexity of the system and the extreme parameters 
involved. However, an observational fact about the depth of the zonal flows of Jupiter and 
Saturn will likely be established given the gravity and magnetic experiments from the Juno 
mission [Bolton, 2010] and the upcoming Cassini Grand Finale [Spilker et al, 2014]. In this 
paper, we focus on the gravity held. 

The physical principle of the gravitational sounding of giant planet zonal flows is not 
complicated; zonal flows will induce local and non-local density perturbations, as well as 
global shape change of the planet, all of which will contribute to perturbations to the external 
non-spherical gravity held. Apart from observational issues such as data coverage, the anal¬ 
ysis of the actual gravity measurement is further complicated by the fact that the background 
external non-spherical gravity held caused by the background uniform rotation is not known 
a priori. Even if one would like to analyze the problem in real space (e.g. directly assess the 
gravity held g rather than the gravity moments 7„), one is still forced to analyze the truncated 
gravity held associated with high-degree gravity moments only (e.g. Ag associated with J„, 
n > 12 for Jupiter and Saturn), in order to retrieve information about the contribution of dif¬ 
ferential rotation to the gravity held. Since the gravity measurements at Jupiter and Saturn 
will not be sensitive to an infinite series of high-degree gravity moments due to the geomet¬ 
ric decay, the accuracy of the individual high-degree gravity moments associated with zonal 
flows from a forward model is an important issue. 

The thermal wind equation (TWE) under the anelastic approximation can be used to 
calculate the gradient of local density perturbations Vp' associated with zonal flows, when 
the zonal flows are much slower than the background rotation. The measured differential ro¬ 
tation on the surface of Jupiter and Saturn is small compare to the background planetary ro¬ 
tation. In terms of the Rossby number Ro - ufQ.oRp (u is the velocity measured in the coro¬ 
tating frame, Qq is the background rotation rate, Rp is the planetary radius), Ro at Jupiter is 
smaller than 0.01, and Ro at Saturn is smaller than 0.05. However, the applicability of TWE 
to further calculate the individual gravity moments for an oblate planet is not guaranteed a 
priori, given the non-local nature of the gravity moments. This applicability has been ac¬ 
tively debated in the recent literature [cf. Kaspi et al, 2010; Kong et al, 2014; Zhang et al, 
2015; Kaspi et al, 2016; Galanti et al, 2017]. 
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In this paper, we present a critical evaluation of the applicability of the thermal wind 
equation to calculate the zonal flow gravity moments. The gravity moments associated with 
deep zonal flows calculated from two versions of the thermal wind equation are compared to 
the full solution to the Euler equation obtained from the Bessel method [Hubbard, 1975; 
Wisdom, 1996; Hubbard, 1999] and the Concentric Maclaurin Spheroid (CMS) method 
[Hubbard, 2013]. A suite of barotropic wind profiles are evaluated. We found that wind in¬ 
duced gravity moments calculated from TWE are in overall agreement with the full solution 
to the Euler equation. However, only when the non-sphericity of the background density and 
effective gravity is taken into account, is the individual high-degree gravity moment calcu¬ 
lated from the thermal wind equation a good approximation to the full solution. Our analysis 
thus suggests that, when analyzing the zonal flow gravity moments of Jupiter and Saturn us¬ 
ing the thermal wind equation, the non-sphericity of the background state should be retained. 

The paper is organized as following, section 2 introduces the definition and properties 
of gravity moments, section 3 presents a detailed comparison of the Euler equation and the 
thermal wind equation, section 4 presents the gravity moments of a uniformly rotating planet 
with polytrope of index unity calculated from the Euler equation using the Bessel method 
and the CMS method, section 5 presents the gravity moments of a differentially rotating 
planet with polytrope of index unity calculated from the Euler equation, the thermal wind 
equation with spherical background state, and the thermal wind equation with non-spherical 
background state, section 6 presents an analysis of what the thermal wind equation misses, 
section 7 summarizes the results and discusses the implications for analyzing the gravity 
measurements of the Juno mission and the Cassini Grand Finale. 

2 Definition and Properties of the Gravity Moments 

The axisymmetric gravity moments J„ are determined by the planetary interior density 
distribution through 



( 1 ) 


in which M is the mass of the planet, a is a reference radius usually chosen to be the mea¬ 
sured equatorial radius of the planet, r is the spherical-radial distance to the center of mass of 
the planet, are the Legendre polynomials of degree n, 0 is the co-latitude measured from 
the spin-axis and the integration is over the entire volume of the planet. 
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It should be immediately realized that 1) if the density distribution is spherically sym¬ 
metric, dp!do - 0, all J„ with n > 1 would be zero; 2) if the density distribution is equatori- 
ally symmetric, p(r, 0) - p(r, n - 0), all odd-degree would be zero. 

If mass, equatorial radius, rotation rate and the gravity moments are the only measure¬ 
ments we have about a planetary body, the interpretation of gravity moments depend on extra 
assumptions/knowledge and forward modeling of the density distributions inside the planet. 
The extra assumptions/knowledge include the composition, temperature, and the equation of 
state (EOS) of the relevant material. The appropriate forward model for the density distribu¬ 
tion inside a fluid planet, for a given EOS, is nothing but the appropriate governing equations 
of fluid dynamics. In the inviscid limit, this set of governing equation is the Euler equation. 
Even though the dynamics being considered are usually simple (e.g. uniform rotation or dif¬ 
ferential rotation on cylinders only), the Euler equations for this particular application are not 
easy to solve due to the fact that we are dealing with self-gravity. The problem is non-local: 
gravity at any local position depends on the density distribution over the entire planet. Math¬ 
ematically, one needs to deal with integro-differetial equations in general. 


3 From the Euler Equation to the Thermal Wind Equation 
3.1 The Euler Equation 


The structure and dynamics of a self-gravitating fluid body in steady-state must be in 
force-balance. This force-balance in the inviscid limit is described by the steady-state Euler 
equation. In an inertial frame, the Euler equation reads 


VP 

(u ■ v)u =-vy„, 

P 


( 2 ) 


in which u is the velocity in the inertial frame, P is the pressure, p is the density, and Vg is 
the gravitational potential, the negative gradient of which is the gravitational acceleration: 




(3) 


Considering self-gravity only, the gravitational potential is determined by the global 
density distribution through 

in which G is the gravitational constant, and the integration is over the entire domain of the 
planet. 
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Under the barotropic assumption (density depends on pressure only) and a velocity 
field that does not violate the barotropic assumption, the density distribution can be deter¬ 
mined entirely from the Euler equation (2) given the total mass and the specific equation of 
state. If the fluid is baroclinic, an additional equation governing the evolution of temperature 
or entropy is needed to determine the density distribution. 

3.2 The Euler Equation in an Inertial Erame for a Uniformly Rotating Planet 

For a uniformly rotating planet, the velocity in the inertial frame reads 

uo = sin 6^, (5) 

here Qq is the constant angular velocity and s is the cylindrical radial distance from the spin 
axis {s - r sin0). 

It can be easily shown that 

(uo ■ V)uo = VQo (6) 

in which Qo is the familiar centrifugal potential 

Qo^-nls'ds'^ -1-. (7) 

The Euler equation now reads 

VP = -pV(Ug + Qo) = -pVU, (8) 

where U is the effective potential defined as U - Vg + Qo- 

When coupled with a specific equation of state (EOS), the solution of the above Euler 
equation yields the shape and internal density distribution of a uniformly rotating planet. It is 
appropriate to denote the properties satisfying equation (8) as the background properties with 
a subscript 0, so equation (8) now reads 

VPo = -poV(yg„ + Qo) = -poVt/o. (9) 


3.3 The Euler Equation for a Planet with Differential Rotation and the Thermal 
Wind Equation 

Now consider a planet with differential rotation, the velocity in the inertial frame reads 

ui = [Qo + Q'(s, z))s^ = QiS(^, (10) 
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here Q' is the angular velocity of the differential rotation, Qi is the total angular velocity 
measured in the inertial frame, while 


u' = ui - uo = z)s^ 


(11) 


is the zonal velocity measured in the non-inertial frame rotating at angular velocity Qq. 


If the angular velocity of the differential rotation depends only on the cylindrical radius 
(do.' Idz - 0), it can be shown that 

(ui ■ V)ui = Va (12) 


here Q is the generalized centrifugal potential 


Q^- f ^\s'ds' = - /" [Qo + n'(s')]^s''^s'- 

Jo Jo 


(13) 


For an arbitrary flow, it can be shown that, 

(ui ■ V)ui = VQo + 2Qoz X u' + (u' ■ V)u', (14) 

it should be recognized that -VQo and -2flo2 x u' are simply the centrifugal acceleration 
and the Coriolis acceleration in the rotating frame with angular velocity Qq. 

Substitute equation (14) and equation (9) into the Euler equation (2), and write the den¬ 
sity and pressure as the sum of the background and the perturbation 

pi^po+p' (15) 


Pi^Po + P\ 


we get 


2poQoz X u' -I- po(u' ■ V)u' -H poVVg' = 
-p'Vy^o-p'VQo-VF' 

- 2p'Qoz X u' - p'(u' ■ V)u' - p'VVgs 


where Vg' is the gravitational potential associated with the density perturbations p' 


y,Kr) 




3 |r - r' 


-p\r')cf‘ 


(16) 

(17a) 

(17b) 

(17c) 


( 18 ) 


An order of magnitude analysis yields the first estimate about the relative importance 
of each of the terms in equation (17). 1) The ratio between the advection term associated 


-7- 



Confidential manuscript submitted to JGR-PIanet.' 

with the zonal flows po(u' ■ V)u' and the Coriolis term 2poQ^z x u' is on the order of the 
Rossby number: ~ 1% for Jupiter and ~ 5% for Saturn. 2) The gravity anomaly associated 
with the density perturbation, VVg/, is smaller than the background gravity by a factor of 
p'!p. Zhang et al. [2015] pointed out that poWg/ could be comparable to p'^Vg^. 3) The 
gradient of the pressure perturbation VP' is likely comparable to p'VUq. This can be shown 
easily for a polytropic equation of state P = The perturbative pressure can now be 

expressed as a function of the density 

P'=P-Po (19a) 

(19b) 

Assuming p'/po ^ 1, one can Taylor expand the above equation. Retaining the first order 
term only, we get 

P' - (1 + -)Kp]l’^p'. (20) 

n “ 

Taking the gradient of P', and making use of the hydrostatic balance of the background state, 
we get 

VP' = -Vvt/o + (l + -)Jifpy"Vp'. (21) 

n n ^ 

(Note, the second term in the RHS of (21) can be comparable to the first term ^p'Vt/o given 
the likely small characteristic scale of p'.) And 4) the ratio of each of the last three terms 
in equation (17) to its corresponding LHS term is |p'/po|. One cautionary note about the 
value of |p'/po|. As we will see, although it is true that |p'/po| is smaller than one in the 
bulk interior of the planet, this is not true for regions very near the surface. 

Taking the curl of the equation (17), retaining only the first term on the left hand side 
(LHS) and the first two terms on the right hand side (RHS), and making use of the mass con¬ 
tinuity equation under the anelastic approximation, we arrive at the generic thermal wind 
equation (TWE) 

(2^0 ■ V)(pou) = -Vp' X geff, (22) 

where the background effective gravity ge® is 

geff = -Vt/o = -V(yg0 + go). (23) 


Note, first, that since the curl has been taken, information has been lost. Specifically, 
there are solutions to the curl free part of equation which can contribute to density perturba¬ 
tions. In addition, there will be density perturbations with non-local origin because of 1) the 
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gravity resulting from the local density anomalies that are required by the TWE and 2) the 
global shape change associated with the net angular moment of the zonal flows. The non¬ 
local density perturbations associated with the gravity anomaly resulting from the local den¬ 
sity perturbations required by the TWE was recognized by Zhang et al. [2015]. However, as 
we will show in section 6, this type of non-local density perturbation contributes very little to 
the high-degree gravity moments. 

One further simplification to the generic thermal wind equation (22) usually adopted 
in estimating the zonal flow gravity field is to assume that the background effective gravity is 
spherically symmetric [e.g. Kaspi et al, 2010; Liu et al, 2013]. One argument for this sim¬ 
plification is the uniqueness of calculated under this assumption despite a non-uniqueness 
in the density perturbations calculated from the TWE. However, we will show that the math¬ 
ematical uniqueness gained from this assumption is not worth the physical relevance being 
sacrificed. And the mathematical non-uniqueness in the density perturbations from the TWE 
can be treated through physically reasonable assumptions. 


4 Gravity Moments of a Uniformly Rotating Planet with Polytrope of Index Unity 

To compare the density perturbations and gravity moments associated with the deep 
zonal flows calculated from the thermal wind equation to the full solution of the Euler equa¬ 
tion, we first solve the Euler equation for a uniformly rotating planet. Here, we adopt a poly¬ 
tropic equation of state 

P = (24) 


in which K is a constant, and the polytropic index n is set to 1. A polytrope of index unity 
not only is a reasonable approximation to the adiabatic equation of state under Jupiter con¬ 
ditions but also makes the Euler equation easier to deal with since VP/p now reduces to 
2KVp. The divergence of the Euler equation (9) yields 


9 2ttG 

V Po + ~^Po = ^ (■ 


AK ^ 


(25) 


which governs the background density distribution subject to the boundary condition that the 
outer boundary defined by pair, 0) -Qis also an equipotential surface 


Uoipoir, 6) - 0) - const. 


(26) 
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As shown by Hubbard [1975], Wisdom [1996], and Hubbard [1999], the general solu¬ 
tion to the above equation takes the form 

^ ma x 

P = Pp+ ^ A„jnikr)P„icos6), (27) 

n=0 

where pp - Q^/lnG, j„ are the spherical Bessel functions of the first of kind of degree n, 
k - ^InGjK, and P„ are the Legendre polynomials of degree n, and A„ are the coefficients 
to be determined by the boundary conditions as well as the total mass. 

The existence of the analytical form of the general solution enables a non-perturbative 
approach to this problem [Hubbard, 1975; Wisdom, 1996; Hubbard, 1999]. Under this cir¬ 
cumstance, there is no need to define the level surfaces and solve for the figure equations ex¬ 
plicitly. Here we point out a few key aspects of this non-perturbative approach which we call 
the Bessel method following Wisdom [1996] and Wisdom and Hubbard [2016]: 1) the coef¬ 
ficients An are the only variables that need to be solved explicitly, both the outer boundary 
shape, which defines the solution domain, and the internal density distribution are uniquely 
determined by A„ ; 2) the outer boundary is not constrained to be an exact ellipsoid of revolu¬ 
tion, and the resulting outer boundary indeed differs from an exact ellipsoid of revolution; 3) 
the traditional geophysical expansion of the external gravitational potential Uq is used to cal¬ 
culate the potential at the outer boundary, given that its convergence under Jupiter and Saturn 
like surface distortions has been shown by Hubbard et al. [2014]. 

The Concentric Maclaurin Spheroid (CMS) method by Hubbard [2013] is a non-perturbative 
method capable of solving for the equilibrium internal density distribution of a rotating planet 
with an arbitrary equation of state. Interested readers should refer to Hubbard [2013] for 
the details of the method. Here we summarize a few key aspects of the CMS method: 1) the 
modeled planet is discretized into a finite number of concentric constant-density spheroids; 

2) the shape of each constant-density spheroid is found iteratively via requiring it to be an 
equipotential surface; 3) the gravitational potential at each level surface is the sum of the 
contributions from every spheroid; 4) the density of each constant-density spheroid needs to 
be adjusted iteratively to match the prescribed equation of state and the fixed total mass of 
the planet. 

Fig. 1 shows the gravity moments and the internal density distribution from these 
two methods with K - 2x10^ [m^kg~^s~^]. The total mass and the background rota¬ 
tion period have been fixed to 1.8983 x Kfi''kg and 9.925 hours respectively, very close to 
the measured value of Jupiter. Tab. 1 compares solutions to the Euler equation (the equa- 
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Table 1. Gravity Moments of a Uniformly Rotating Polytrope of Index Unity Planet from Different Meth¬ 
ods. 




Polytrope Model^ 

Jupiter" Bessel Method 

Polytrope Model^ 

CMS Method (521) 

Fractional 

Difference 


Mass [kg] 

1.89819 X 10^^ 

1.8983 X lO^’^ 

1.8983 X 1022 


Rotation Period [hrs] 

9.925 

9.925 

9.925 


Equatorial Radius [km] 

71492 

71418.75 

71407.15 

1.62 X 10-4 

Polar Radius [km] 

66854 

66875.28 

66868.54 

1.01 X 10-4 

J2 X 10“® 

14696.43 

13948.95 

13953.32 

3.13 X 10-4 

74 X 10“® 

-587.14 

-528.81 

-529.14 

6.18 X 10-4 

76 X 10“® 

34.25 

29.86 

29.89 

9.15 X 10-4 

78 


-2.108 X 10“® 

-2.110 X 10-° 

1.20 X 10-2 

7io 


1.716 X 10“’^ 

1.718 X 10-2 

1.48 X 10-2 

7i2 


-1.541 X 10“^ 

-1.544 X 10-* 

1.76 X 10-2 

7i4 


1.488 X 10“'^ 

1.491 X 10-° 

2.02 X 10-2 

7i6 


-1.517 X 10-1° 

-1.520 X 10-1° 

2.28 X 10-2 

7i8 


1.614 X 10““ 

1.618 X 10-11 

2.54 X 10-2 

720 


-1.777 X 10-12 

-1.782 X 10-12 

2.78 X 10-2 


“ The physical values are from Archinal et al. [2009] & Jacobson [2003] with G = 6.67408 x 10 
^ P = Kp^ with = 2 X 10^ [m^kg~^s~^]. 


torial radius, the polar radius, and the gravity moments up to degree 20) from the Bessel 
method and the CMS method for a uniformly rotating planet with a polytrope of index unity. 
It can be seen from Tab. 1 that the fractional differences of the gravity moments, defined as 
|y„(Bessel)-7„(CMS)|/|7„(Bessel)|, are on the order of 3x 10“^ ~ 3x10“^. This agrees with 
the Wisdom and Hubbard [2016] assessment that the discretization error for a 512-spheroid 
CMS is on the order of 1 x 10“'^ ~ 1 x 10“^. The observed values of Jupiter are listed as well. 
It can be seen that the first three gravity moments of this model planet are reasonably close to 
those measured at Jupiter. 

The right panel of Fig. 2 shows the effective ellipticity of the equipotential surface as a 
function of the equatorial radius of the equipotential surface. It can be seen that the effective 
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Spehrical Harmonic Degree 



Figure 1. Gravity moments J„ and internal density distribution of a uniformly rotating planet with fixed 
mass, rotation rate, and a polytropic equation of state of index unity with K = 2 x 10^ [m^kg~^ s~^]. For J„, 
filled (open) symbols represent positive (negative) values. For J„, solutions to the full Euler equation from the 
Bessel method and the Concentric Maclaurin Spheroid method (with 521 spheroids) are shown here. Due to 
the excellent agreement between the two methods, the two solutions appear indistinguishable on this plot. 


ellipticity decreases from ~ 0.35 near the surface of the planet to ~ 0.29 near the center of 
the planet, and this shape change occurs mostly in the outer part of the planet. The effective 
ellipticity is defined as e(r) - - r^lr^, where and are the equatorial radius and the 

polar radius of level surfaces. We call this effective ellipticity due to the fact that the level 
surface are not exact ellipsoids of revolution. The left panel of Fig. 2 shows the deviation of 
the outer boundary equipotential surface and mid-shell equipotential surface from an exact 
ellipsoid of revolution with the equatorial radius and polar radius fixed to the corresponding 
values of the equipotential surface. The deviation is dominated by sin^ 28 with an amplitude 
~ 5 X 10“^ at the outer boundary, and thus corresponds to the second order correction in 
the standard expansion of level surface in terms of the effective ellipticity (e.g. equation 30.3 
in Zharkov and Trubitsyn 1978). We noticed that some published solutions of this problem 
are based on the assumption that the outer boundary shape is an exact ellipsoid of revolution 
[e.g. Kong et al, 2013, 2015]. 

5 Gravity Moments of a Differentially Rotating Planet with Polytrope of Index 
Unity 

We now turn to a planet with differential rotation. The total mass, background rota¬ 
tion rate, and the equation of the state are taken to be the same as those in the study of a uni- 
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Deviation of Equipotential Surfaces 
from an Exact Ellipsoid of Revolution 


Internal Shape 




Equatorial Radius of Equipotential Surface [ap] 


Figure 2. Deviation of equipotential surfaces from that of an exact ellipsoid of revolution with the equato¬ 
rial radius and the polar radius fixed to the corresponding values of the equipotential surface, and the internal 
shape of equipotential surfaces measured by the effective ellipticity e=^\ — r^lr^. The deviation of the 
outer boundary from an exact ellipsoid of revolution shows a dominant component of sin^ 29 with amplitude 
~ 5 X 10“^. This corresponds to the second order correction in the standard expansion of level surface in 

terms of the effective eccentricity (e.g. equation 30.3 in Zharkov and Trubitsyn 1978). 


formly rotating planet. Note that the equatorial radius and the shape of a differentially ro¬ 
tating planet would be different from those of a uniformly rotating planet. The differential 
rotation is chosen to have angular velocity as a function of cylindrical-radial distance only 
(dQ.ldz = 0). This consideration is mainly motivated by the fact that such a velocity pro¬ 
file does not violate the barotropic assumption, and a full solution to the Euler equation can 
be obtained without further complications. (Of course, there is no reason to suppose that the 
planets obey this precisely.) The full solution of the Euler equation can then be compared to 
that obtained from the thermal wind equation. 

Fig. 3 shows a series of zonal wind profiles adopted in this study. All the zonal wind 
profiles feature a prominant band of equatorial super rotation. Different wind profiles are 
characterized by the different half-amplitude-width (HAWD), which is defined as the cylin¬ 
drical radius at which the amplitude of the zonal wind equals half of the peak amplitude. To 
ensure convergence, the generalized centrifugal potential (13) is approximated by a polyno¬ 
mial expansion of cylindrical radius s and truncated at degree 24 following Hubbard [1982]; 
Kaspi et al. [2016]; Wisdom and Hubbard [2016]; Galanti et al. [2017]. The density per¬ 
turbations and gravity moments associated with this wind profile are then calculated using 
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Distance to the Spin-Axis [ap] 



Figure 3. Zonal wind profiles considered in this study. All the zonal winds are assumed to be constant 
along the direction parallel to the spin-axis. Different wind profiles are characterized by the different half- 
amplitude-width (HAWD), which is defined as the fractional cylindrical radius at which the amplitude of the 
zonal wind equals half of the peak amplitude. 


four different approaches: the Bessel method (full solution to the Euler equation), the CMS 
method (full solution to the Euler equation), the thermal wind equation with spherical back¬ 
ground state, and the thermal wind equation with non-spherical background state. 


5.1 Euler Equation Solution from the Bessel Method and the CMS Method 

With differential rotation on cylinders and a polytrope of index unity, the divergence of 
the Euler equation (2) now reads 


o 2nG 9 Q 


(28) 


where Q is the generalized centrifugal potential (13). 


Eor the Bessel method, the solution to this equation takes the same functional form as 
(27). The difference is that the specific solution pp is now a function of cylindrical radius s 
rather than a constant. The specific solution pp(s) [e.g. Hubbard, 1999] for the wind profile 
we are considering can be obtained via numerically integrating equation (28) with the inner 
boundary condition pp(0) = Qg/27rG. This inner boundary condition is only valid when the 
coefficient of the term in the polynomial expansion of Q = is set to 1 /2 Qq. 

After obtaining ppis), we can solve for the internal density distribution and gravity moments 
Jn associated with the differential rotation using the same non-perturbative approach. 
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Table 2. Zonal wind induced AJ„ calculated from Bessel method and CMS method for HAWD = 0.80. 



[10-6] 

Bessel Method CMS (521) 

1 Bessel-CMS (521) | 

Fractional Difference 


A 72 

46.045 

46.290 

2.45 X 10“' 

3.30 X 10- 

AJa 

-16.861 

-16.929 

6.71 X 10-2 

7.22 X 10- 

AJe 

4.856 

4.873 

1.62 X 10-2 

1.25 X 10- 

AJs 

-0.510 

-0.508 

2.65 X 10-2 

4.52 X 10- 


A7io 

-0.394 

-0.398 

4.35 X 10-2 

1.84 X 10-2 

A7i2 

-0.203 

-0.204 

1.01 X 10-2 

5.24 X 10-2 

A7i4 

0.0202 

0.0211 

8.67 X 10-4 

4 X 10-2 

A7i6 

-0.0575 

-0.0582 

6.41 X 10-4 

1.11 X 10-2 

A7i8 

0.0146 

0.0146 

6.35 X 10-2 

4.33 X 10-2 

A72o 

0.0105 

0.0108 

2.75 X 10-4 

2.61 X 10-2 


For the CMS method, one only needs to replace the centrifugal potential term Qo in the 
equipotential surface equations with the generalized centrifugal potential term Q. There is no 
need to obtain a specific solution pp in the CMS method. 

The left panel of Fig. 4 compares the total gravity moments associated with the 
zonal flows shown in Fig. 3 to the gravity moments associated with the background uniform 
rotation, while the right panel of Fig. 4 compares the wind-induced gravity moments only, 
defined as A7„ = 7„ (Uniform Rotation + Wind) - 7„ (Uniform Rotation). The solution ob¬ 
tained from the Bessel method and the CMS method are in very good agreement and appear 
almost identical in the figure, thus only the solution from the Bessel method is shown here 
for clarity. Tables 2-4 presented the tabulated comparison between the Bessel method and 
the CMS method for the wind induced gravity moments. It can be seen from the tables that 
these two methods agree very well for all three different wind profiles; 1) the absolute differ¬ 
ences in AJn are smaller than 3 x 10“® beyond degree 10, and 2) the fractional differences in 
Jr, are on the order of 3%, except for Jn with absolute values smaller than 3 x 10“^. In con¬ 
trast, Kaspi et al. [2016] reported much larger discrepancies between the CMS method and 
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Table 3. Zonal wind induced AJ„ calculated from Bessel method and CMS method for HAWD = 0.90. 



[10-6] 

Bessel Method CMS (521) 

1 Bessel-CMS (521) | 

Fractional Difference 


A 72 

5.835 

5.882 

4.69 X 10-2 

3.17 X 10-"^ 

AJa 

-4.662 

-4.679 

1.75 X 10-2 

6.46 X 10-"^ 

AJe 

2.495 

2.512 

1.74 X 10-2 

1.38 X 10-2 

AJs 

-1.147 

-1.157 

1.04 X 10-2 

3.96 X 10-2 


A7io 

-0.397 

-0.400 

3.40 X 10-2 

6.42 X 10-2 

A7i2 

-0.0439 

-0.0431 

7.50 X 10-4 

1.22 X 10-2 

A7i4 

-0.0695 

-0.0716 

2.06 X 10-2 

3.03 X 10-2 

A7i6 

0.0670 

0.0686 

1.65 X 10-2 

2.47 X 10-2 

A7i8 

-0.0281 

-0.0288 

7.40 X 10-4 

2.64 X 10-2 

A72o 

-0.00237 

-0.00224 

1.36 X 10-4 

5.74 X 10-2 


the Bessel method. In Kaspi et al. [2016], the absolute differences in high-degree J„ are typ¬ 
ically on the order of 1 x 10“^, and the fractional differences in high-degree are typically 
on the order of 30%. We are confident in the good agreement we obtained between the CMS 
method and the Bessel method, given that Wisdom and Hubbard [2016] have independently 
shown the very good agreement between these two methods for the uniform rotation case. 

It should be emphasized here that in both the Bessel method and the CMS method, we 
are solving the full Euler equation to get the total gravity moments for the case with zonal 
flows, rather than solving for the perturbations AJ„ only. This is fundamentally different 
from the approach taken by Kong et al. [2014, 2015], in which the wind induced gravity is 
treated as a perturbation. One aspect of the full solution is that the outer boundary shape of 
the planet gets further changed when zonal winds are included. 

It can be seen from the left panel of Fig. 4 that the gravity moments associated with 
the zonal flows shown in Fig. 3 exceed those from the background rotation by more than 
100% starting around degree 10. For the low-degree gravitational moments, the wind in¬ 
duced contribution increases as the depth of the wind increases (see Fig. 4 & Tab. 2 - 4). For 
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Table 4. Zonal wind induced AJ„ calculated from Bessel method and CMS method for HAWD = 0.975. 



[10-6] 

Bessel Method CMS (521) 

1 Bessel-CMS (521) | 

Fractional Difference 


A 72 

1.963 

1.929 

3.40 X 10-2 

3.11 X 10-"^ 

A74 

-0.380 

-0.343 

3.77 X 10-2 

5.47 X 10-"^ 

A76 

0.189 

0.169 

1.69 X 10-2 

3.84 X 10-"^ 

A78 

-0.128 

-0.121 

7.15 X 10-2 

2.06 X 10-2 


A7io 

0.0929 

0.0904 

3.40 X 10-2 

8.45 X 10-2 

A7i2 

-0.0685 

-0.0683 

2.45 X 10-4 

2.60 X 10-2 

A7i4 

0.0488 

0.0495 

7.09 X 10-4 

1.42 X 10-2 

A7i6 

-0.0374 

-0.0383 

9.14 X 10-2 

2.44 X 10-2 

A7i8 

0.0263 

0.0270 

6.18 X 10-4 

2.35 X 10-2 

A72o 

-0.00733 

-0.00754 

2.11 X 10-4 

2.87 X 10-2 


the high-degree gravitational moments, however, the wind induced contributions are sim¬ 
ilar for the three different wind depth considered here (see Fig. 4 & Tab. 2 - 4). This fea¬ 
ture should caution about the attempt to estimate the wind-induced perturbations to the low- 
degree gravity moments based on future measurements of the high-degree gravity moments. 

5.2 Thermal Wind Equation with Spherical Background Density and Gravity 

We now proceed to solve for the gravity moments and density perturbations associ¬ 
ated with the same zonal flows using the thermal wind equation. We first consider the ther¬ 
mal wind equation under the simplification that reduces both the background density and 
the background effective gravity to a spherically symmetric state. This is the simplification 
adopted in many published calculations of the zonal flow gravity using the thermal wind 
equation [e.g. Kaspi et al., 2010; Liu et al, 2013]. The thermal wind equation now reads 


2|Gol 


g[po(r)u'] 

dz 


1 dp' 
~r~^ 


ee X (-|go(r)|er). 


(29) 
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Gravity Moments from Euler Equation 




Figure 4. Gravity moments associated with deep equatorial zonal flows calculated from the Euler equation. 
Jn are shown in the left panel, while Ay„ = (Differential Rotation) - (Uniform Rotation) are shown in the 
right panel. For and Ay„, filled (open) circles represent positive (negative) values. 



Figure 5. Gravity moments associated with deep equatorial zonal flows calculated from the Euler equation 
compared to those calculated from the spherical thermal wind equation and the non-spherical thermal wind 
equation for HAWD=0.90. Only Ai„ are shown here, and filled (open) circles represent positive (negative) 
values. 


The spherically symmetric background density and background effective gravity correspond¬ 
ing to a planet with the same mass and the same polytropic index unity equation of state can 
be obtained analytically. 
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Figure 6. Gravity moments associated with deep equatorial zonal flows calculated from the Euler equation 
compared to those calculated from the spherical thermal wind equation and the non-spherical thermal wind 
equation for HAWD=0.80 & HAWD=0.975. Only /S.J„ are shown here, and filled (open) circles represent 
positive (negative) values. 


One motivation for adopting a spherically symmetric background state in the thermal 
wind equation to calculate the wind induced gravity moments is the uniqueness of 7,, despite 
the non-uniqueness of wind induced density perturbations. With spherically symmetric back¬ 
ground state, the TWE yields the gradient of the density perturbations Vp' along the 9 direc¬ 
tion dp'! do. To get the density perturbations p'(r, 9) with spherically symmetric background 
state, one would integrate dp' jdO along the 9 direction 

p'(r,61)= / —rd9' + p'^{r), (30) 

7e'=o on 

here p'^ir) is a “constant of integration" which is a function of r only. This "constant of in¬ 
tegration" resulting from spherical thermal wind equation makes zero contribution to the 
gravity moments with n > 1, since 



p'^{r)r”Pn{cos9)(f’r - 0, 


n > 1. 


(31) 


The gravity moments associated with the zonal flows shown in Fig. 3 calculated from 
the thermal wind equation with spherically symmetric background density and spherically 
symmetric background gravity is compared to those calculated from the Euler equation in 
Fig. 5 & 6 and in Tab. 5. Only AJ„ are shown in Fig. 5 & 6 and in Tab. 5. It can be seen 
that 1) overall, the solution from the spherical TWE is in order-of-magnitude agreement with 
the full solution; 2) however, some of the individual high-degree gravity moments calculated 
from this simplified thermal wind equation can be wrong by more than 100% and can take 
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Table 5. Zonal wind induced AJ„ for HAWD=0.90 calculated from four methods. 


[10-6] 

Bessel 

CMS 

Spherical 

Non-spherical 


Method 

Method 

TWE 

TWE 

AJi 

5.835 

5.882 

5.836 

5.897 

AJ4 

-4.662 

-4.679 

-3.586 

-3.826 

AJe 

2.495 

2.512 

1.810 

2.098 

A/g 

-1.147 

-1.157 

-0.728 

-0.978 

A7io 

0.397 

0.400 

0.155 

0.333 

A/i2 

-0.0439 

-0.0431 

0.0682 

-0.0268 

A7i4 

-0.0695 

-0.0716 

-0.0945 

-0.0681 

A7i6 

0.0670 

0.0686 

0.0486 

0.0614 

A7i8 

-0.0281 

-0.0288 

-0.0001 

-0.0242 

AJiq 

-0.00237 

-0.00224 

-0.01786 

-0.00338 


the wrong sign. This is consistently the case for the series of zonal wind profiles adopted in 
this study. 

5.3 Thermal Wind Equation with Non-spherical Background Density and Effec¬ 
tive Gravity 

We now proceed to calculate the density perturbation and gravity moments associated 
with zonal flows using the generic thermal wind equation with non-spherical background 
density distribution and non-spherical background effective gravity 

(2^0 ■ V)(poU) = -Vp' X geff. (32) 


For a polytrope of index of unity, the background effective gravity geu is simply 

geff = -Vt/o = —- - 2KVpo, (33) 

Po 

which can be easily calculated since Vpo is entirely determined by the coefficients A„. 
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With non-spherical effective gravity, the TWE now yields the gradient of the density 
perturbations Vp' along the tangent of equipotential surfaces in the meridional plane instead 
of the gradient of the density perturbations along the 6 direction. To get the density perturba¬ 
tions p'(r, 6) with non-spherical effective gravity, one would need to integrate Vp' along the 
tangent of equipotential surfaces in the meridional plane 



(34) 


here ^ is measured along the direction perpendicular to the equal-potential surface, I and dl' 
are the meridional arc length measured on the equal-potential surface, and p'^(^) is a “con¬ 
stant of integration" which is a function of ^ instead of r. This "constant of integration" re¬ 
sulting from non-spherical thermal wind equation makes a non-zero contribution to the grav¬ 
ity moments. Since Vp^(^) x geff = 0, this “constant of integration" has zero contribution 
to the generic thermal wind equation. It is clear then that TWE itself cannot supply the "con¬ 
stant of integration". However, since p'^(^) is only a function of the gravity moments as¬ 
sociated with this “constant of integration" should be proportional to the background 7„, and 
the pre-factor should be on the order of the ratio of the wind-induced mass anomaly to the 
total mass of the planet which is a very small quantity. Since the wind-induced 7„ are orders 
of magnitude larger than the background J„ beyond n - 10, the correction to by p^Cf) 
must be negligible beyond n = 10. Thus we set p'^(^) to zero in all our non-spherical TWE 
calculations. 

It can be seen from Eig. 5 & 6 and Tab. 5 that for n > 10, A7„ associated with zonal 
flows calculated from the thermal wind equation with non-spherical background state is 
much closer to the full solution than those calculated from the thermal wind equation with 
spherical background state. All A7„ greater than 3 x 10“® now take the correct sign, and the 
amplitude difference is now within 50% for individual A7„. 

6 What the Thermal Wind Equation Misses: Global Shape Change, Non-local and 
Irrotational Density Contrihutions 

As discussed in section 3.3, density perturbations that contribute irrotationally to the 
Euler equation and that have a non-local origin must necessarily exist. Eig. 7 compares the 
density perturbations associated with deep zonal flows (the case with HAWD=0.90 in Eig. 

3) calculated from the Euler equation and the thermal wind equation. The density perturba¬ 
tions calculated from two version of the thermal equation are both localized and are visu¬ 
ally similar, thus only the density perturbations calculated from the thermal wind equation 
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Density Perturbation in the Equatorial Plane 




Figure 7. Density perturbations associated with deep zonal flows (the case with HAWD=0.90 in Fig. 3) 
calculated from the Euler equation and the thermal wind equation. The leftmost panel shows the density 
perturbations in the equatorial plane while the right two panels show the meridional cut. It can be seen that 
the thermal wind equation captures the local wind-induced density perturbations but misses the non-local 
large-scale density perturbations. 


with spherical background state is shown for clarity. It can be seen from Fig. 7 that the ther¬ 
mal wind equation captures the local density perturbations directly associated with the local 
zonal flows but misses the large-scale density perturbations (with a dominant degree-2 struc¬ 
ture) associated with the global shape change of the planet related to the net positive angular 
moment of the zonal flows. The zonal wind profiles we are considering all have net positive 
total angular moment. As a result the planet shrinks in the polar direction while it expands in 
the equatorial direction. 

Fig. 8 shows the outer boundary shape change of the entire planet caused by the local¬ 
ized zonal flows shown in Fig. 3. It can be seen that the equatorial radius increased by ~ 25 
km, while the polar radius decreased by ~ 4 km for an equatorial super rotation with half am¬ 
plitude width at 0.80 /?/. Small-scale shape changes spatially correlated with the local zonal 
flows are also evident Fig. 8. 

Zhang et al. [2015] proposed a gravitational correction to the thermal wind equation 
by including the gravitational anomalies from the local density perturbation required by the 
TWE; the curl of the po^Vg' term. We performed an independent calculation of the gravi¬ 
tational anomaly correction as proposed in Zhang et al. [2015], and find that for the high- 
degree gravity moments beyond n = 12, the gravitational anomaly term can only provide 
corrections on the order of 1%. This is consistent with the independent analysis of Galanti 
etal. [2017]. 
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Global Outer Boundary Shape Change 
Caused by Localized Zonal Winds 



Figure 8. Global outer boundary shape change induced by the localized deep zonal flow shown in Fig. 3. It 
can be seen that the localized equatorial jet with half amplitude width (HAWD) of 0.80 can induce a change 
to the outer boundary position by ~ 25 km near the equator and by ~ 4 km near the poles where there are no 
local zonal flows. 


A cautionary note about the wind-induced low-degree gravity moments; they are in¬ 
significant corrections to the background low-degree gravity moments (e.g 0.3% correction 
to J 2 and 3% to ^4 for the deepest wind profile considered here with HAWD = 0.80/?/; 
0.01% correction to J 2 and 0.07% correction to Ja, for HAWD = 0.915Rj), and can be easily 
offset by uncertainties in the background state (such as uncertainties in our knowledege about 
the equation of state, thermal state, heavy element distribution, etc.) 

7 Summary and Discussion 

In this paper, we present a critical examination of the applicability of the thermal wind 
equation under anelastic assumption to calculate the gravity moments associated with deep 
zonal flows of giant planets. We first derive the thermal wind equation from the Euler equa¬ 
tion and show that the thermal wind equation is a good approximation to the local dynam¬ 
ics when the Rossby number of the zonal flow measured in the co-rotation frame is much 
smaller than one. It is also pointed out that the thermal wind equation is a local treatment 
even when the background effective gravity is used, while the full problem is non-local. 

We then solve the full Euler equation for a rotating self-gravitating planet with poly¬ 
trope of index unity. The Bessel method [Hubbard, 1975; Wisdom, 1996; Hubbard, 1999] 
and the Concentric Maclaurin Spheroid method [Hubbard, 2013] are employed. We first 
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solve for the shape, density distribution, and the gravity moments of a uniformly rotating 
planet. It is shown that the outer boundary shape has a significant deviation, on the second 
order, from an exact ellipsoid of revolution. The impact of the assumption that the outer 
boundary shape is an exact ellipsoid of revolution adopted in some studies of this problem 
[e.g. Kong et al, 2013, 2015] on the solutions of gravity moments requires further investiga¬ 
tion. 


For Jupiter-like zonal flows but confined to the equatorial region and assumed to be 
constant on cylinders (e.g. Fig. 3), the associated density perturbations and gravity moments 
are then calculated from four different methods: the Bessel method, the CMS method, the 
thermal wind equation with spherical background state, and the thermal wind equation with 
non-spherical background state. The full solutions to the Euler equation obtained from the 
Bessel method and the CMS method show excellent agreement. 

Concerning the applicability of the thermal wind equation, we And that 1) overall, the 
solution from the spherical TWE is in order-of-magnitude agreement with the full solution; 
2 ) a few individual high-degree gravity moments calculated from the spherical thermal wind 
equation can be wrong by 100% and can take the wrong sign; 3) the individual high-degree 
gravity moments calculated from the thermal wind equation with non-spherical background 
density and non-spherical effective gravity are a good approximation to the full solution, the 
difference is within 50%; 4) for low-degree gravity moments associated with zonal flows, 
global shape change to the planet caused by the net angular moments of the zonal flows is 
important. This global shape change is missed in the thermal wind equation as well as in the 
thermal-gravitational wind equation [Zhang ef a/., 2015]. However, the wind-induced low- 
degree gravity moments may not be a concern since they are most likely indiscernible from 
uncertainties in the background state. 

Eor baroclinic zonal winds with velocity variations along the direction of spin-axis, 
we don’t know how to solve the full Euler equation to obtain the density perturbation, shape 
change, and gravity moments. Eor baroclinic winds, the equipotential surface and equal- 
density surface are mis-aligned in general. Eurthermore, the generalized centrifugal poten¬ 
tial Q (e.g. equations 12 - 13) cannot be defined for a baroclinic flow. To be more specific, 

(u ■ V)u of a baroclinic flow cannot be written as a gradient of a scalar potential. The current 
technique to solve the Euler equation, such as the Bessel method and the CMS method, itera¬ 
tively And the solution via requiring the equal-density surface being an equipotential surface. 
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thus cannot be straightforwardly generalized to deal with baroclinic flows. Non-perturbative 
methods are yet to be developed to solve the full equation for baroclinic zonal flows. The 
thermal wind equation, on the other hand, remains a valid approximation for low-Rossby 
flows. Furthermore, the thermal wind equation, with spherical or non-spherical background 
state, requires significantly less computational resources than the full Euler equation. 

A clear message from this study to the analysis of the gravity measurements from Juno 
and Cassini is that to calculate the wind-induced gravity moments using the thermal wind 
equation, taking the non-spherical nature of the background density and effective gravity into 
account would yield much more accurate representation of the full solution. 
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